# -*- coding: utf-8 -*-
""" A Python interface module to some GDAL library functions to translate
    LAS Variable Record content (Geotiff keys) into WKT or Proj4 strings

    :Author:
      - 20120816-20120929 Roberto Vidmar

    :Revision:  $Revision: 8 $
                $Date: 2013-02-01 10:12:57 +0100 (Fri, 01 Feb 2013) $

    :Copyright: 2012-2012
                Nicola Creati <ncreati@inogs.it>
                Roberto Vidmar <rvidmar@inogs.it>

    :License: MIT/X11 License (see :download:`LICENSE.txt
                               <../../LICENSE.txt>`)
"""

import ctypes
import sys
import struct
from lasExceptions import ConversionFailed, UnimplementedException

SHORT = 'h'
DOUBLE = 'd'

# Tag types (from geotiff.h)
TAGTYPE_BYTE     = 1
TAGTYPE_SHORT    = 2
TAGTYPE_LONG     = 3
TAGTYPE_RATIONAL = 4
TAGTYPE_ASCII    = 5
TAGTYPE_FLOAT    = 6
TAGTYPE_DOUBLE   = 7
TAGTYPE_SBYTE    = 8
TAGTYPE_SSHORT   = 9
TAGTYPE_SLONG    = 10
TAGTYPE_UNKNOWN  = 11

pinfo_t = ctypes.c_ushort

#===============================================================================
class GTIFDefn(ctypes.Structure):
  MAX_GTIF_PROJPARMS = 10
  _fields_ = [
    ('Model',             ctypes.c_short),
    ('PCS',               ctypes.c_short),
    ('GCS',               ctypes.c_short),
    ('UOMLength',         ctypes.c_short),
    ('UOMLengthInMeters', ctypes.c_double),
    ('UOMAngle',          ctypes.c_short),
    ('UOMAngleInDegrees', ctypes.c_double),
    ('Datum',             ctypes.c_short),
    ('PM',                ctypes.c_short),
    ('PMLongToGreenwich', ctypes.c_double),
    ('Ellipsoid',         ctypes.c_short),
    ('SemiMajor',         ctypes.c_double),
    ('SemiMinor',         ctypes.c_double),
    ('TOWGS84Count',      ctypes.c_short),
    ('TOWGS84',           ctypes.c_double * 7),
    ('ProjCode',          ctypes.c_short),
    ('Projection',        ctypes.c_short),
    ('CTProjection',      ctypes.c_short),
    ('nParms',            ctypes.c_int),
    ('ProjParm',          ctypes.c_double * MAX_GTIF_PROJPARMS),
    ('ProjParmId',        ctypes.c_int * MAX_GTIF_PROJPARMS),
    ('MapSys',            ctypes.c_int),
    ('Zone',              ctypes.c_int),
    ('DefnSet',           ctypes.c_int),
    ]

  def __repr__(self):
    """ A representation of all fields of a :class:`GTIFDefn` instance
    """
    # I love strings!
    s = "%s:\n" % self.__class__.__name__
    for f in self._fields_:
      val = self.__getattribute__(f[0])
      p = str(val).partition('Array')[-1]
      if p:
        size = p.split()[0].strip('_')
        s += "  %s = %s\n" % (f[0], repr([val[i] for i in xrange(int(size))]))
      else:
        s += "  %s = %s\n" % (f[0], val)
    return s

#===============================================================================
class TIFFMethod(ctypes.Structure):
  _fields_ = [
    ('get', ctypes.c_void_p),
    ('set', ctypes.c_void_p),
    ('type', ctypes.c_void_p)]

#===============================================================================
class GeoKey(ctypes.Structure):
  _fields_ = [
    ('gk_key', ctypes.c_int),
    ('gk_size', ctypes.c_size_t),
    ('gk_type', ctypes.c_int),
    ('gk_count', ctypes.c_long),
    #('gk_data', ctypes.c_char_p),
    ('gk_data', ctypes.POINTER(ctypes.c_char)),
    ]

#===============================================================================
class GTiff(ctypes.Structure):
  """ The class to interface Python code with libgdal. We need to fill this
      structure to get back later a WKT string
  """
  GeoKeyDirectoryTag = 34735
  GeoDoubleParamsTag = 34736
  GeoAsciiParamsTag = 34737

  _fields_ = [
    ('gt_tif', ctypes.c_void_p),
    ('gt_methods', TIFFMethod),
    ('gt_flags', ctypes.c_int),
    ('gt_version', pinfo_t),
    ('gt_rev_major', pinfo_t),
    ('gt_rev_minor', pinfo_t),
    ('gt_num_keys', ctypes.c_int),
    ('gt_keys', ctypes.POINTER(GeoKey)),
    ('gt_keyindex', ctypes.POINTER(ctypes.c_int)),
    ('gt_keymin', ctypes.c_int),
    ('gt_keymax', ctypes.c_int),
    ('gt_short', ctypes.POINTER(pinfo_t)),
    ('gt_double', ctypes.POINTER(ctypes.c_double)),
    ('gt_nshorts', ctypes.c_int),
    ('gt_ndoubles', ctypes.c_int)]

  def __init__(self, geokeys, geoAsciiParams=None, geoDoubleParams=None):
    """ Create a new instance of the structure and populate it with data
        from Geotiff like "geokeys" tuple, geoAsciiParams string and
        geoDoubleParams tuple from LAS Variable Length Records.

        This is an example::

          [
          (1, 1, 0, 13),
          (1024, 0, 1, 1),
          (2048, 0, 1, 32767),
          (2050, 0, 1, 32767),
          (2056, 0, 1, 7030),
          (3072, 0, 1, 32767),
          (3074, 0, 1, 32767),
          (3075, 0, 1, 1),
          (3076, 0, 1, 9001),
          (3080, 34736, 1, 1),
          (3081, 34736, 1, 0),
          (3082, 34736, 1, 3),
          (3083, 34736, 1, 4),
          (3092, 34736, 1, 2)
          ]
          "UTM Zone 15, Northern Hemisphere|NAD83|"
          (0.0, -81.0, 0.9996, 500000.0, 0.0)
  """
    #for gk in geokeys:
    #  print gk
    #print geoDoubleParams

    # Create an empty list of c_doubles, we will need the REFERENCE to them
    # Do *NOT* delete the next statement
    self.c_doubles = []
    key0 = geokeys[0]
    keys = [k[0] for k in geokeys][1:]
    self.gt_methods = TIFFMethod()
    self.gt_short = (pinfo_t * (len(geokeys) * 4))()
    self.gt_version = key0[0]
    self.gt_rev_major = key0[1]
    self.gt_rev_minor = key0[2]
    self.gt_num_keys = key0[3]
    self.gt_keys = (GeoKey * len(geokeys))()
    self.gt_keymin = min(keys)
    self.gt_keymax = max(keys)
    self.gt_nshorts = (len(geokeys) * 4) # Every key has four shorts

    # Create an array of 32768 elements to find the geokey tuple
    # from its key: 32767 mean unknown *AND* INDEX 0 means key not found! 
    self.gt_keyindex = (ctypes.c_int * 32768)()

    if geoDoubleParams is None:
      self.gt_double = (ctypes.c_double * 1)()
      self.gt_ndoubles = 0
    else:
      self.gt_double = (ctypes.c_double * len(geoDoubleParams))()
      self.gt_ndoubles = len(geoDoubleParams)

      for i, d in enumerate(geoDoubleParams):
        self.gt_double[i] = ctypes.c_double(d)

    geoKeysFlattened = sum(geokeys, ())
    for i, v in enumerate(geoKeysFlattened):
      self.gt_short[i] = ctypes.c_ushort(v)

    for j, key in enumerate(geokeys[1:]):
      i = j + 1 # AAAAARGH!!!!!
      self.gt_keys[i].gk_key = key[0]
      self.gt_keyindex[key[0]] = i
      if key[1] == 0:
        # Data type is unsigned short
        self.gt_keys[i].gk_size = 2
        self.gt_keys[i].gk_type = TAGTYPE_SHORT
        ptr = key[3]
      elif key[1] == GTiff.GeoDoubleParamsTag:
        # Data type is double
        self.gt_keys[i].gk_size = 8
        self.gt_keys[i].gk_type = TAGTYPE_DOUBLE
        # 3 days of HARD WORK for these three lines
        d = geoDoubleParams[key[3]]
        self.c_doubles.append(d)
        value = ctypes.c_double(d)
        ptr = ctypes.pointer(value)
      elif key[1] == GTiff.GeoAsciiParamsTag:
        # Data type is ascii
        self.gt_keys[i].gk_size = 1
        self.gt_keys[i].gk_type = TAGTYPE_ASCII
        # Remove null character
        string = geoAsciiParams[key[3]: key[3] + key[2] - 1]
        ptr = ctypes.create_string_buffer(string)
      else:
        raise UnimplementedException(
          "wTIFFTagLocation value '%d' is not implemented." % key[1])

      self.gt_keys[i].gk_count = key[2]
      self.gt_keys[i].gk_data = ctypes.cast(ptr, ctypes.POINTER(ctypes.c_char))

  def __repr__(self):
    """ Unambiguous representation of the GeoKeys, GeoDoubleParams and
        geoAsciiParams
    """
    s = "%s:\n" % self.__class__.__name__
    geoKeys, geoDoubles, geoAscii = self.geoKeys()
    s += "Geokeys: (\n"
    for gk in geoKeys:
      s += "  %s\n" % repr(gk)
    s += ")\ngeoDoubleParams:\n"
    s += "  %s\n" % repr(geoDoubles)
    s += "geoAsciiParams:\n"
    s += "  %s\n" % repr(geoAscii)
    return s

#-------------------------------------------------------------------------------
# Public methods
#-------------------------------------------------------------------------------

  def geoKeys(self):
    """ Return (GeoKeys, GeoDoubleParams, GeoAsciiParams) tuple

    :returns: (GeoKeys, GeoDoubleParams, GeoAsciiParams)
    :rtype: tuple
    """
    keys = [self.gt_keys[i].gk_key for i in xrange(1, self.gt_num_keys)]
    # Keys must be sorted first
    keys.sort()

    # Build the geoKeys
    geoDoubleParams = ()
    geoAsciiParams = ''
    geoKeys = ((self.gt_version, self.gt_rev_major,
      self.gt_rev_minor, self.gt_num_keys - 1), )

    for k in keys:
      index = self.gt_keyindex[k]
      key = self.gt_keys[index]
      if key.gk_type == TAGTYPE_DOUBLE:
        double = struct.unpack(DOUBLE, key.gk_data[:8])[0]
        geoDoubleParams += (double, )
        geoKeys += ((
          key.gk_key,
          GTiff.GeoDoubleParamsTag,
          key.gk_count,
          len(geoDoubleParams) - 1), )
      elif key.gk_type == TAGTYPE_SHORT:
        # ANOTHER 3 DAYS OF CRAZY WORK TO UNDERSTAND THIS!
        # Look at gdal-1.9.0/frmts/gtiff/libgeotiff/geo_get.c
        s = ctypes.cast(
          ctypes.pointer(key.gk_data), ctypes.POINTER(ctypes.c_char))
        short = struct.unpack(SHORT, s[:2])[0]
        geoKeys += ((
          key.gk_key,
          0,
          key.gk_count,
          short), )
      elif key.gk_type == TAGTYPE_ASCII:
        ascData = key.gk_data[:key.gk_count]
        geoKeys += ((
          key.gk_key,
          GTiff.GeoAsciiParamsTag,
          key.gk_count,
          len(geoAsciiParams)), )
        geoAsciiParams += ascData

    return geoKeys, geoDoubleParams, geoAsciiParams

#===============================================================================
class FILE(ctypes.Structure):
  pass

pFile = ctypes.POINTER(FILE)

#===============================================================================
# Link SOME GDAL Library functions
#lib = ctypes.CDLL('C:\Python27\Lib\site-packages\osgeo\gdal19.dll')
#lib = ctypes.CDLL('/usr/local/lib/libgdal.so.1')
libpath = ctypes.util.find_library('gdal')
lib = ctypes.CDLL(libpath)
#lib = ctypes.CDLL('/d0/Downloads/gdal-1.9.0-Bobo/.libs/libgdal.so.1')

# Link GTIFGetDefn
GTIFGetDefn = lib.GTIFGetDefn
GTIFGetDefn.argtypes = [ctypes.POINTER(GTiff), ctypes.POINTER(GTIFDefn)]
GTIFGetDefn.restype = ctypes.c_int

# Link GTIFPrintDefn
GTIFPrintDefn = lib.GTIFPrintDefn
GTIFPrintDefn.argtypes = [ctypes.POINTER(GTIFDefn), pFile]
GTIFPrintDefn.restype = ctypes.c_int

# Link GTIFGetOGISDefn
GTIFGetOGISDefn = lib.GTIFGetOGISDefn
GTIFGetOGISDefn.argtypes = [ctypes.POINTER(GTiff), ctypes.POINTER(GTIFDefn)]
GTIFGetOGISDefn.restype = ctypes.c_char_p

# Link GTIFSetFromOGISDefn
GTIFSetFromOGISDefn = lib.GTIFSetFromOGISDefn
GTIFSetFromOGISDefn.argtypes = [ctypes.POINTER(GTiff),
  ctypes.POINTER(ctypes.c_char)]
GTIFSetFromOGISDefn.restype = ctypes.c_int

# Link GTIFNew
GTIFNew = lib.GTIFNew
GTIFNew.argtypes = [ctypes.c_void_p]
GTIFNew.restype = ctypes.POINTER(GTiff)

#-------------------------------------------------------------------------------
# Link `Python File` as `File`
PyFile_AsFile = ctypes.pythonapi.PyFile_AsFile
PyFile_AsFile.argtypes = [ctypes.py_object]
PyFile_AsFile.restype = pFile

#-------------------------------------------------------------------------------
# Get stdout
stdout_file = PyFile_AsFile(sys.stdout)

#-------------------------------------------------------------------------------
def geoKeysFromWkt(wkt):
  """ Return geoKeys tuple, geoDoubleParams and geoAsciiParams from wkt string

  :returns: (GeoKeys, GeoDoubleParams, GeoAsciiParams)
  :rtype: tuple
  """
  # Create e new (empty) GTiff structure
  pNewGtif = GTIFNew(ctypes.c_void_p())

  # Set wkt string to the new object
  rc = GTIFSetFromOGISDefn(pNewGtif, wkt)
  if not rc:
    raise ConversionFailed(
      "GTIFSetFromOGISDefn failed with WKT = '%s'" % wkt)

  gtif = pNewGtif.contents
  return gtif.geoKeys()

#-------------------------------------------------------------------------------
def normalize(geoKeyDirectoryTag, geoDoubleParamsTag=None,
  geoASCIIParamsTag=None):
  """ Return "normalized" (i.e. translated into WKT) version of geotiff tags
      found in variable length records of a LAS file

      See gdal-1.9.0/frmts/gtiff/libgeotiff/geo_normalize.c

  :param geoKeyDirectoryTag: Mandatory GeoKeyDirectoryTag
  :type geoKeyDirectoryTag: :class:`VLRecord` instance
  :param geoDoubleParamsTag: GeoDoubleParamsTag
  :type geoDoubleParamsTag: :class:`VLRecord` instance
  :param geoASCIIParamsTag: GeoASCIIParamsTag
  :type geoASCIIParamsTag: :class:`VLRecord` instance
  :returns: wkt
  :rtype: string
  """
  geoKeys = geoKeyDirectoryTag.getData()
  if geoDoubleParamsTag:
    geoDoubleParams = geoDoubleParamsTag.getData()
  else:
    geoDoubleParams = None
  if geoASCIIParamsTag:
    geoAsciiParams = geoASCIIParamsTag.getData()
  else:
    geoAsciiParams = None

  # Create GTiff tags
  gt = GTiff(geoKeys, geoAsciiParams=geoAsciiParams,
    geoDoubleParams=geoDoubleParams)

  # Normalize GTiff tags
  defn = GTIFDefn()
  rc = GTIFGetDefn(ctypes.byref(gt), ctypes.byref(defn))

  #print defn
  if not rc:
    raise ConversionFailed(
      "GTIFGetDefn failed converting '%s'" % repr(gt))

  # Get wkt string from gtif structure (GTiff tags)
  wkt = GTIFGetOGISDefn(ctypes.byref(gt), ctypes.byref(defn))
  return wkt

#-------------------------------------------------------------------------------
if __name__ == '__main__':
  # The Geo Keys:
  geoKeys = [
  (1, 1, 0, 13),       # a Version 1 GeoTIFF GeoKey, Rev. 1.0, 13 Keys
  (1024, 0, 1, 1),     # TIFFTagLocation = 0 count = 1 Projected = 1
  (2048, 0, 1, 32767),
  (2050, 0, 1, 32767),
  (2056, 0, 1, 7030),
  (3072, 0, 1, 32767),
  (3074, 0, 1, 32767),
  (3075, 0, 1, 1),
  (3076, 0, 1, 9001),
  (3080, 34736, 1, 1), # GeoDoubleParamsTag = 34736, count = 1 index = 1
  (3081, 34736, 1, 0), # GeoDoubleParamsTag = 34736, count = 1 index = 0
  (3082, 34736, 1, 3), # GeoDoubleParamsTag = 34736, count = 1 index = 3
  (3083, 34736, 1, 4),
  (3092, 34736, 1, 2)
  ]
  geoDoubles = (0.0, -81.0, 0.9996, 500000.0, 0.0)
  gtif = GTiff(geoKeys, geoDoubleParams=geoDoubles)
  # Normalize GTiff tags
  defn = GTIFDefn()
  rc = GTIFGetDefn(ctypes.byref(gtif), ctypes.byref(defn))
  if not rc:
    raise ConversionFailed(
      "GTIFGetDefn failed converting '%s'" % repr(gtif))

  # Get wkt string from gtif structure (GTiff tags)
  wkt = GTIFGetOGISDefn(ctypes.byref(gtif), ctypes.byref(defn))
  print "\n\nWKT from geokeys and GeoDoubles:\n", wkt
  print "Creating now a new object from this WKT string...",

  # Create e new (empty) GTiff structure
  pNewGtif = GTIFNew(ctypes.c_void_p())

  # Set wkt string to the new object
  rc = GTIFSetFromOGISDefn(pNewGtif, wkt)
  if not rc:
    raise ConversionFailed(
      "GTIFSetFromOGISDefn failed with WKT = '%s'" % wkt)

  # Normalize GTiff tags
  defnNew = GTIFDefn()
  rc = GTIFGetDefn(pNewGtif, ctypes.byref(defnNew))
  if not rc:
    raise ConversionFailed(
      "GTIFGetDefn failed converting '%s'" % repr(defnNew))

  # It worked !
  print "done."
  print "\nThis is the new object structure printed by GTIFPrintDefn:"
  GTIFPrintDefn(ctypes.byref(defnNew), stdout_file)
  wktNew = GTIFGetOGISDefn(pNewGtif, ctypes.byref(defnNew))
  print "\nwkt == wktNew is '%s'" % str(wkt==wktNew)

  geoKeysFromWkt(wktNew)
